capture log close

set more off
set matsize 8000 
program drop _all
set scheme s1mono

log using "JMBAL_JCR_Replication.log", replace
	cd ""		/*Add file path for your machine*/
	
	
*PURPOSE: 	This file will replicate the models and graphs reported in Braithwaite
*			and Licht, "The Effect of Civil Society Organizations and Democratization
*			Aid on Civil War Onset". Before running the file:
*				1) Save JMBAL_JCR_Replication.dta to your machine
*				2) Change the directory command <cd> at the top of the do-file to 
*				the folder where you have saved JMBAL_JCR_Replication.dta
*				3) Delete first line of code below instructions if outreg2 is installed
*NOTE: 		The simulation (beginning line 380) requires extensive computing power.  
*			It may take a very long time to run. On my machine (quad core 3.3Ghz proc, 
*			solid state, 24G RAM) processing time for the simulations is 5 hours. If 
*			preferred, users may comment out the simulations and use 
*			JMBAL_SimulationResults.dta to make Figure 2.
*AUTHOR:	Amanda A. Licht, Binghamton University
*MACHINE:	LICHTBINGOFFICE, Lenovo Thinkstation running StataMP 13
*DATE:		10/23/2019
*MODIFIED:		
*FILE:		JMBAL_JCR_Replication.do
*LOG:		Logs/JMBAL_JCR_Replication.log
*DATA:		JMBAL_JCR_Replication.dta, 
*Sources:	Tierney, Michael J., Daniel L. Nielson, Darren G. Hawkins, J. 
*			Timmons Roberts, Michael G. Findley, Ryan M. Powers, Bradley Parks, 
*			Sven E. Wilson, and Robert L. Hicks. 2011. More Dollars than Sense: 
*			Refining Our Knowledge of Development Finance Using AidData. 
*			World Development 39 (11): 1891-1906. 
*			Gleditsch, Kristian S. 2002. "Expanded Trade and GDP Data," 
*			Journal of Conflict Resolution 46: 712-24.
*			Coppedge, Michael, John Gerring, Staffan I. Lindberg, Svend-Erik 
*			Skaaning, Jan Teorell, David Altman, Michael Bernhard, M. Steven Fish, 
*			Adam Glynn, Allen Hicken, Carl Henrik Knutsen, Kyle Marquardt,
*			Kelly McMann, Farhad Miri, Pamela Paxton, Daniel Pemstein, Jeffrey Staton, 
*			Eitan Tzelgov, Yi-ting Wang, and Brigitte Zimmerman. 2016. �V-Dem 
*			[Country-Year/Country-Date] Dataset v6.2.� Varieties of Democracy (V-Dem) Project.
*			Cheibub, Jos� Antonio, Jennifer Gandhi, and James Raymond Vreeland. 2010. 
*			�Democracy and Dictatorship Revisited.� Public Choice, vol. 143, no. 2-1, pp. 67-101.
*			Melander, Erik and Ther�se Pettersson & Lotta Themn�r (2016) Organized violence, 
*			1989-2015. Journal of Peace Research 53(5).
*			Gleditsch, Nils Petter, Peter Wallensteen, Mikael Eriksson, Margareta Sollenberg, 
*			and H�vard Strand (2002) Armed Conflict 1946-2001: A New Dataset. Journal of Peace 
*			Research 39(5).
*			World Bank,  2017, "Historical Classifications by Income", 
*			available online at <https://datahelpdesk.worldbank.org/knowledgebase/
*			articles/906519-world-bank-country-and-lending-groups>. Accessed 3/7/17.
*			Bailey, Michael and Strezhnev, Anton and Voeten, Erik, (2013).
*			"Estimating Dynamic State Preferences from United Nations Voting
*			Data".Available at SSRN: http://ssrn.com/abstract=2330913
*			(September 25, 2013).
*OUTPUT:	
*	Tables:	JMBAL Double Endog 2sls.txt
*			JMBAL Double Endog IV Probit.txt
*	Data: 	JMBAL_RobustAnalysisIV.dta		- includes lincom of 2sls combined coefs	
*			JMBAL_ObservedValuesSims		- Observed values of control variables for
*												simulation procedure
*			JMBAL_InterpIVprob				- 1K sims of mean predicted probability and
*												marginal effects over observed values of
*												covariates from double endog IVprobit 
*			JMBAL_InterpIVprobCIs	- CIs around above using Normality assumption
*	Graphs:	JMBAL_FIGURE1.gph
*			JMBAL_FIGURE2.gph

*Install outreg2 for making tables
net install outreg2

use  JMBAL_JCR_Replication.dta, clear


xtset ccode year

*********************************************************************	
*1: Run 2sls to get fit stats and instrument assessments
*********************************************************************

#delimit ;
ivregress 2sls	onset2v414 	L1.v2x_cspart
								L1.lfoodaid_rusd 
								L1.lbudgetaid_rusd
								Llpop
								Grgdppc
								peacetime*
				(L1.lgovtcivsoc_rusd partyXdemaid_rusd= L1.UNGAdist_USA L1UNGAdist_USAXL1v2xcspart L1.UNGAdist_mu L1UNGAdist_muXL1v2x_cspart) 
				if aidelg==1, 
				robust cluster(ccode)  first;														
#delimit cr

*Save parameter matrices
mat def b2sls = e(b)
mat def V2sls = e(V)

*get instrument tests:
estat firststage, forcenonrobust

*get combined coefficients for table:
sum L1.v2x_cspart
lincom L1.lgovtcivsoc_rusd  + partyXdemaid_rusd*`r(mean)'
	local ccAidIV = round(`r(estimate)', .001)
	local ccAIDIVse = round(`r(se)', .001)
	
sum L1.lgovtcivsoc_rusd 
lincom L1.v2x_cspart + partyXdemaid_rusd*`r(mean)'
	local ccParty = round(`r(estimate)', .001)
	local ccPartyse = round(`r(se)', .001)
	
	*Write to table:
	#delimit ;
	outreg2 using "Double Endog 2sls rusd.txt",
		dec(3) ctitle(Interactive Model)
		addtext("Combined Coefficient Aid, `ccAidIV'(`ccAIDIVse')" 
				"Combined Coefficient Civil Society, `ccParty'(`ccPartyse')") 
		replace;	
	#delimit cr		

	
*********************************************************
*Quick interpretation of interaction effect w/lincom
*********************************************************
*make variables to hold range of conditioning variable values:
sum  L1.lgovtcivsoc_rusd
	local min =`r(min)'
	local max = `r(max)'
gen AidX =	`r(min)' in 1
	replace AidX = AidX[_n-1] + ((`max'-`min')/19) if _n>1&_n<=20
		sum AidX

sort AidX		
	sum  L1v2x_cspart
	local min =`r(min)'
	local max = `r(max)'
gen csoX =	`r(min)' in 1
	replace csoX = csoX[_n-1] + ((`max'-`min')/19) if _n>1&_n<=20
		sum csoX
		
*make places to store combined coefs:
gen CC_Aid		= .
gen CC_Aid_se 	= .
gen CC_CSO		= .
gen CC_CSO_se 	= .		
	
	
		*Loop lincom commands over 20 values of conditioning variables:
		forvalues num = 1(1)20	{
	
			sort AidX
			local AID = AidX in `num'		/*local value for correct row*/
			local cso = csoX in `num'		
			
			*Calculate combined coef and std error using Delta method
			lincom L1.v2x_cspart + partyXdemaid_rusd*`AID'	
				*Store values in the CC variable in correct row
				replace CC_CSO 		= `r(estimate)' in `num'
				replace CC_CSO_se 	= `r(se)' in `num'
			
			lincom L1.lgovtcivsoc_rusd  + partyXdemaid_rusd*`cso'
				replace CC_Aid		= `r(estimate)' in `num'
				replace CC_Aid_se 	= `r(se)' in `num'				
								}
								
	*Get CIs:
	foreach var in CC_Aid CC_CSO	{
		gen `var'_lo90 = `var' - 1.65*`var'_se
		gen `var'_lo95 = `var' - 1.96*`var'_se
		gen `var'_hi95 = `var' + 1.96*`var'_se
		gen `var'_hi90 = `var' + 1.65*`var'_se
									}
									
	
	
****************************************************************
*FIGURE 1: effects on probability of conflict onset
****************************************************************

*Store Y-axis values for rug plots
 gen rugCSO = -.3
 gen rugAID = -.005
 gen pipe = "|"
 
*Marginal Effect of Civil Society Strength: 
#delimit ;
twoway (line CC_CSO CC_CSO_l* CC_CSO_h* AidX, sort 
			lw(medium thin thin thin thin)
			lp(solid dash dot dot dash)
			lc(black gs10 gs10 gs10 gs10))
		(scatter rugCSO lgovtcivsoc_rusd, sort msym(none) mlab(pipe) mlabpos(0)) 	,
		title("Civil Society Strength", size(medium))
		xtitle("Democracy Aid", size(small))
		yline(0, lc(gs12))
legend(off)		
name(CC_CSO, replace);

*Marginal Effect of Democracy Assistance
#delimit ;
twoway (line CC_Aid CC_Aid_l* CC_Aid_h* csoX, sort 
			lw(medium thin thin thin thin)
			lp(solid dash dot dot dash)
			lc(black gs10 gs10 gs10 gs10))
		(scatter rugAID v2x_cspart, sort msym(none) mlab(pipe) mlabpos(0))	,
		title("Democracy Assistance", size(medium))
		xtitle("Civil Society Strength", size(small))
				yline(0, lc(gs12))
legend(off)
name(CC_AID, replace);	

*Combine the two panels to make one graph:
#delimit ;
graph combine CC_CSO CC_AID, rows(1)
	title("Conditional Effects on Conflict Onset", size(medium))
saving("JMBAL_FIGURE1.gph", replace)	;
#delimit cr 

**********************************************************************
*Additional Instrument Tests
**********************************************************************

*ivregress command doesn't provide F-tests, so I also run the model using ivreg2:
xtset ccode year		
#delimit ;
ivreg2  	onset2v414 	L1.v2x_cspart
								L1.lfoodaid_rusd 
								L1.lbudgetaid_rusd
								Llpop
								Grgdppc
								peacetime*
				(L1.lgovtcivsoc_rusd partyXdemaid_rusd= L1.UNGAdist_USA L1UNGAdist_USAXL1v2xcspart L1.UNGAdist_mu L1UNGAdist_muXL1v2x_cspart) 
				if aidelg==1, 
				robust cluster(ccode)  first;														
#delimit cr



*Save data w/Combined coefs:
save JMBAL_RobustAnalysisIV.dta, replace


************************************************************************
*2: IVPROBIT
************************************************************************
xtset ccode year
#delimit ;
ivprobit  	onset2v414 	L1.v2x_cspart
								L1.lfoodaid_rusd 
								L1.lbudgetaid_rusd
								Llpop
								Grgdppc
								peacetime*
				(L1.lgovtcivsoc_rusd partyXdemaid_rusd= L1.UNGAdist_USA L1UNGAdist_USAXL1v2xcspart L1.UNGAdist_mu L1UNGAdist_muXL1v2x_cspart) 
				if aidelg==1, 
				vce(cluster ccode)  first;														
#delimit cr

*Get combined coefficients for table:
sum L1.v2x_cspart
lincom L1.lgovtcivsoc_rusd  + partyXdemaid_rusd*`r(mean)'
	local ccAidIV = round(`r(estimate)', .001)
	local ccAIDIVse = round(`r(se)', .001)
	
sum L1.lgovtcivsoc_rusd 
lincom L1.v2x_cspart + partyXdemaid_rusd*`r(mean)'
	local ccParty = round(`r(estimate)', .001)
	local ccPartyse = round(`r(se)', .001)

	*Write to table:
	#delimit ;
	outreg2 using "Double Endog IV Probit rusd.txt",
		dec(3) ctitle(Interactive Model)
		addtext("Combined Coefficient Aid, `ccAidIV'(`ccAIDIVse')" 
				"Combined Coefficient Civil Society, `ccParty'(`ccPartyse')") 
		replace;	
	#delimit cr	
	
*Save parameter matrices	
mat def bIV = e(b)
mat def VIV = e(V)


****************************************
*INTERPRETATION
****************************************

*Below, a simulation procedure to generate predicted probabilities across the 
*joint range of CSO strength and democracy aid will be implemented to get the 
*values used in Figure 2.


*First, we need to get dataset of values of control variables:
keep if e(sample)
#delimit;
keep lgovtcivsoc_rusd
	lfoodaid_rusd 
	lbudgetaid_rusd
	Llpop
	Grgdppc
	peacetime;
#delimit cr	

sum lgovtcivsoc_rusd
sum lgovtcivsoc_rusd, detail

*get values of aid to use in sims:
_pctile lgovtcivsoc_rusd if lgovtcivsoc_rusd>0 , p(5, 10, 25, 30, 45, 50, 55, 60, 75, 80, 90, 95)
return list

save JMBAL_ObservedValuesSims.dta, replace



*****************************************	
*Program for IV probit interpretation:
*****************************************

program def ivpred
	syntax, aid(real) cso(real)	/*tell Stata to expect real values called aid and cso*/
		
		*Draw a vector of coefficients from the parameter matrices:
		clear
		set obs 1
		
		drawnorm b1-b42, means(bIV) cov(VIV)
		
		*Expand and merge w/observed values of covariates:
		expand 4106
		merge 1:1 _n using JMBAL_ObservedValuesSims, nogen
		
		*Grab local values of aid and cso, to ease XB calculation 
			gen double aid = `aid'
			gen cso = round(`cso', .001)
		
		*Linear index:
		#delimit ;
		gen XB =	b1*aid 		+
					b2*aid*cso 	+
					b3*cso		+
					b4*lfoodaid_rusd +
					b5*lbudgetaid_rusd +
					b6*Llpop +
					b7*Grgdppc +
					b8*peacetime +
					b9*peacetime^2+
					b10*peacetime^3+
					b11;
		#delimit cr
		
		
		*Get combined coefs:
			gen CCaid = b1+b2*cso
			gen CCcso = b3+b2*aid
			
		*Get probabilities and MEs:
			gen ProbOnset 	= normal(XB)
			gen pdfOnset	= normalden(XB)
			gen MEaid 		= pdfOnset*CCaid
			gen MEcso 		= pdfOnset*CCcso
			
		*collapse to in-sample mean across observed values:
		collapse (mean) ProbOnset MEaid MEcso
		
			*reclaim values of aid and cso strength:
			gen double aid = `aid'
			gen cso = round(`cso', .001)
			
		*Save results:	
		append using JMBAL_InterpIVprob.dta
		save JMBAL_InterpIVprob.dta, replace
		
end	/*close program*/

*Make place to store simulation results:
clear
save 	JMBAL_InterpIVprob.dta, replace emptyok

*MAY TAKE MANY HOURS TO COMPLETE:  comment out if preferred
*Run simulations 1000 times at each combination of CSO and Aid values:	
set more off	
forvalues cso = .02(.097).99	{
forvalues aid = 1(1)20	{
	simulate, reps(1000): ivpred, cso(`cso') aid(`aid')
																}
								}		
	
***************************
*Get Confidence Intervals	
***************************

*Open simulation results:
use JMBAL_InterpIVprob.dta,  clear

*Use collapse to get average and standard deviation:
#delimit ;
collapse (mean) ProbOnset MEaid MEcso 
		(sd) ProbOnset_sd=ProbOnset
			MEaid_sd = MEaid
			MEcso_sd = MEcso, by(cso aid);
#delimit cr

*Calculate CIs using Normality assumption:
foreach var in ProbOnset MEaid MEcso	{
	gen `var'_lo90 = `var' - 1.65*`var'_sd
	gen `var'_lo95 = `var' - 1.96*`var'_sd
	gen `var'_hi95 = `var' + 1.96*`var'_sd
	gen `var'_hi90 = `var' + 1.65*`var'_sd
										}


*variables to tag subset of the variables' values:
sort cso aid
browse cso aid 	
bysort cso: gen aidlev = _n 
sort aid cso
bysort aid: gen csolev = _n 

tab1 aidlev csolev
gen csolev2= 1 if csolev==1
	replace csolev2= 2 if csolev==6
	replace csolev2= 3 if csolev==11
	
gen aidlev2 = 1 if aidlev==1
		replace aidlev2 = 2 if aidlev==7
		replace aidlev2 = 3 if aidlev==12
	
	lab def aidfmt2 1 "Minimum" 2 "Median" 3 "95th %tile"
	lab val aidlev2 aidfmt2
	
	lab def csofmt 1 "Minimum" 2 "Average" 3 "Maximum"
	lab val csolev2 csofmt
	
	label var aid "Democracy Assistance"
	label var cso "CSO Participation Index"

save JMBAL_InterpIVprobCIs.dta, replace


***********************************
*FIGURE 2:
***********************************

*Summarize and save min and max values to locals:
sum ProbOnset 
	local min = round(r(min), .01)
	local max = round(r(max), .01)
*contour plot:
#delimit ;
twoway 	(contour ProbOnset cso aid, ccuts(`min'(.01)`max') 
			ecolor(dkgreen)  crule(int)
			ztitle("Predicted Probability of Onset", orient(rvertical))
			ytitle("CSO Participation") ylabel(,angle(horiz) labsi(small))
			xtitle("Democracy Aid", size(small))
			
			title("Predicted Probability of Civil Conflict Onset" 
			"by CSO Strength and Democracy Aid", size(medium))),
name(Pr1, replace)
saving(JMBAL_FIGURE2.gph, replace)	;
#delimit cr

log close
